# clear labels for versions
versions <- tibble(
version = c("v1", "v2", "v3", "v4", "v5", "v6", "v7"),
vlabel = factor(
c( "$Priors\\,Do\\,Not\\,Vary\\,by\\,State\\,and\\,Date$",
"$\\beta$ Centered at Emp. Value",
"$P(S_1|untested)$ Centered at Emp. Value",
"$P(S_1|untested)$ and $\\beta$ Centered at Emp. Values",
"$\\beta$ Centered at Emp. Value (Spline Smoothed)",
"$P(S_1|untested)$ and $\\beta$ Centered at Emp. Values ($\\beta$ Spline Smoothed)",
"$\\,$CTIS Priors that Do Not Vary by State or Date$\\,$"
),
levels = c(
"$Priors\\,Do\\,Not\\,Vary\\,by\\,State\\,and\\,Date$",
"$\\beta$ Centered at Emp. Value",
"$P(S_1|untested)$ Centered at Emp. Value",
"$P(S_1|untested)$ and $\\beta$ Centered at Emp. Values",
"$\\beta$ Centered at Emp. Value (Spline Smoothed)",
"$P(S_1|untested)$ and $\\beta$ Centered at Emp. Values ($\\beta$ Spline Smoothed)",
"$\\,$CTIS Priors that Do Not Vary by State or Date$\\,$"
) )
)
# read state-level results
targets_dir <- here("_targets_6_26")
state_res <- tar_read(state_v1,
store =targets_dir) %>%
bind_rows(
tar_read(state_v2,
store =targets_dir)
) %>%
bind_rows(
tar_read(state_v3,
store =targets_dir)
) %>%
bind_rows(
tar_read(state_v4,
store =targets_dir)
)%>%
bind_rows(
tar_read(state_v5,
store =targets_dir)
)%>%
bind_rows(
tar_read(state_v6,
store =targets_dir)
) %>%
bind_rows(
tar_read(state_v7,
store =targets_dir)
)
covidestim <- tar_read(covidestim_biweekly_state,
store =targets_dir)
dates <- readRDS(here("data/data_raw/date_to_biweek.RDS"))
codes <- read_csv(here('data/demographic/statecodes.csv'))
state_res <- state_res %>%
rename(state=fips) %>%
left_join(versions) %>%
mutate(state=toupper(state)) %>%
left_join(covidestim, relationship= "many-to-many") %>%
left_join(dates, relationship = "many-to-many") %>%
rename(name=state_name) %>%
group_by(biweek) %>%
mutate(mindate=min(date),
maxdate=max(date)) %>%
ungroup()
state_res <- state_res %>%
filter(version %in% c("v2", "v5", "v6", "v7"))
theme_c <- function(...){
# font <- "Helvetica" #assign font family up front
theme_bw() %+replace% #replace elements we want to change
theme(
#text elements
plot.title = element_text( #title
size = 16, #set font size
face = 'bold', #bold typeface
hjust = .5,
vjust = 3),
plot.subtitle = element_text( #subtitle
size = 12,
hjust = .5,
face = 'italic',
vjust = 3), #font size
axis.title = element_text( #axis titles
size = 12), #font size
axis.text = element_text( #axis text
size = 9),
legend.text = element_text(size = 12),
# t, r, b, l
plot.margin = unit(c(1,.5,.5,.5), "cm"),
legend.position = "right",
strip.text.x = element_text(size = 11, face = "bold", color="white"),
strip.text.y = element_text(size = 11, face = "bold", color="white",angle=270),
strip.background = element_rect(fill = "#3E3D3D")
) %+replace%
theme(...)
}
Summarizing Concordance with Covidestim
# corrected %>%
# mutate(in_interval = case_when(
# exp_cases_lb <= infections & infections <= exp_cases_ub ~ "In Interval",
# exp_cases_lb > infections ~ "Below Interval",
# exp_cases_ub < infections ~ "Above Interval"
# )) %>%
# filter(!is.na(in_interval)) %>%
# group_by(in_interval, state, vlabel) %>%
# summarize(n_per_county=n()) %>%
# group_by(vlabel, in_interval) %>%
# summarize(n_per_version = sum(n_per_county)) %>%
# group_by(vlabel) %>%
# mutate(total = sum(n_per_version)) %>%
# ungroup() %>%
# mutate(prop=n_per_version/total)
by_version <- state_res %>%
select(-date) %>%
distinct() %>%
mutate(in_interval = case_when(
exp_cases_lb <= infections & infections <= exp_cases_ub ~ "In Interval",
exp_cases_lb > infections ~ "Below Interval",
exp_cases_ub < infections ~ "Above Interval"
)) %>%
filter(!is.na(in_interval)) %>%
group_by(in_interval, vlabel) %>%
summarize(n=n()) %>%
group_by(vlabel) %>%
mutate(total = sum(n)) %>%
mutate(prop=n/total)
labels <- by_version %>%
arrange(prop) %>%
pull(vlabel) %>%
as.character() %>%
unique()
by_version %>%
mutate(in_interval = factor(in_interval,
levels = c("Above Interval",
"In Interval",
"Below Interval"
))) %>%
ggplot(aes(x= fct_reorder(vlabel,prop,max),
fill = in_interval, y =prop)) +
geom_bar(stat="identity",position="stack") +
theme_c() +
coord_flip()+ scale_fill_manual(values=c("In Interval" = "#79D2AF",
"Above Interval" = "#D2AF79",
"Below Interval"="#7997D2"),
breaks = c("Below Interval",
"In Interval",
"Above Interval")) +
labs(x="",
y= "Proportion",
fill = "Covidestim Median:",
title = "Proportion of Simulation Intervals Where Covidestim Median\nWas Within, Above, or Below the Median",
subtitle = "Including All States") +
theme_c() +
theme_c(legend.position="right",
legend.title = element_text(face="bold", size = 15),
legend.text= element_text( size = 15),
legend.spacing.y = unit(3, 'pt'),
axis.text.y = element_text(size = 17, hjust=1)) +
theme(plot.title = element_text(size = 20),
plot.subtitle = element_text(size=18,
face='italic',
margin=margin(4,0,4,0)),
axis.text.x=element_text(size=12),
axis.title = element_text(size=14)) +
scale_x_discrete(labels = (TeX(labels)))
Simulation Intervals Over Time
all_versions <- state_res %>%
group_by(state) %>%
summarize(state_n=n_distinct(version)) %>%
filter(state_n == max(state_n)) %>%
pull(state)
state_res %>%
filter(state %in% sample(unique(all_versions),5)) %>%
ggplot() +
geom_ribbon(aes(x = date,
ymin = exp_cases_lb,
ymax = exp_cases_ub),
alpha = .8,
show.legend=FALSE,
fill="#596281") +
geom_linerange(aes(xmin = mindate,
xmax=maxdate,
y = positive,
color = 'obs'),
size=.5) +
facet_grid(name~vlabel,
labeller=labeller(vlabel =as_labeller(
TeX, default=label_parsed)),
scales="free_y") +
scale_y_continuous(labels = comma) +
scale_x_date(date_breaks = "3 months",
date_labels = "%b %Y") +
theme_c(
legend.position = "top",
) +
scale_fill_manual(values = pal) +
labs(x="",
y = "Estimated Infections",
title = "Simulation Intervals Over Time",
subtitle = "By State") +
scale_color_manual(values = c('obs' = 'red'), labels = 'Positive Tests',
name = '') +
guides(color = guide_legend(override.aes = list(linewidth =2)))
Compare Versions
subset <- state_res %>%
filter(vlabel == "$\\,$CTIS Priors that Do Not Vary by State or Date$\\,$")
pal <- c("#247C90", "red")
names(pal) <- c(unique(subset$vlabel), "Covidestim")
state_res %>%
filter(state %in% sample(unique(all_versions),5)) %>%
filter(version %in% c("v2", "v5", "v6", "v7")) %>%
ggplot() +
geom_ribbon(aes(x = date,
ymin = exp_cases_lb,
ymax = exp_cases_ub,
fill=vlabel),
alpha = .8,
show.legend=FALSE #,
# fill="#596281"
) +
geom_linerange(aes(xmin = mindate,
xmax=maxdate,
y = positive,
color = 'obs'),
size=.5) +
facet_grid(name~vlabel,
labeller=labeller(vlabel =as_labeller(
TeX, default=label_parsed)),
scales="free_y") +
scale_y_continuous(labels = comma) +
scale_x_date(date_breaks = "3 months",
date_labels = "%b %Y") +
theme_c(
legend.position = "top",
) +
scale_fill_manual(values = pal) +
labs(x="",
y = "Estimated Infections",
title = "Simulation Intervals Over Time",
subtitle = "By State") +
scale_color_manual(values = c('obs' = 'red'), labels = 'Positive Tests',
name = '') +
guides(color = guide_legend(override.aes = list(linewidth =2)))
Compare States
subset %>%
ggplot() +
geom_ribbon(aes(x = date,
ymin = exp_cases_lb,
ymax = exp_cases_ub,
fill = vlabel),
alpha = .8,
show.legend=FALSE) +
geom_ribbon(aes(x = date,
fill = 'Covidestim',
ymin = infections.lo,
ymax = infections.hi),
alpha = .6) +
facet_wrap(~state, ncol=5, scales = "free_y",
labeller= labeller(vlabel =as_labeller(TeX, default=label_parsed))) +
scale_y_continuous(labels = comma) +
scale_x_date(date_breaks = "3 months",
date_labels = "%b %Y") +
theme_c(axis.text.x=element_text(size=6),
axis.text.y = element_text(size=8),
legend.position = "top") +
scale_fill_manual(values =pal,
labels = TeX(names(pal)),
name='',
breaks="Covidestim") +
labs(y = "Estimated Infections",
title = paste0("Estimated Infections by State")) +
guides(color = guide_legend(override.aes = list(linewidth =2),
nrow=2,
byrow=TRUE))
ggsave(here('figures/intervals-and-covidestim-all-states.pdf'))
Ratio of Estimated to Observed
subset %>%
ggplot() +
geom_ribbon(aes(x = date,
ymin = exp_cases_lb/positive,
ymax = exp_cases_ub/positive,
fill = vlabel),
alpha = .9,
show.legend=FALSE,
fill= "#247C90") +
facet_wrap(~state, ncol=5,
labeller= labeller(vlabel =as_labeller(TeX, default=label_parsed))) +
scale_y_continuous(labels = comma) +
scale_x_date(date_breaks = "3 months",
date_labels = "%b %Y") +
theme_c(axis.text.x=element_text(size=6),
axis.text.y = element_text(size=8),
legend.position = "top") +
scale_fill_manual(values =pal,
labels = TeX(names(pal)),
name='',
breaks="Covidestim") +
labs(y = "Ratio of Estimated Infections to Observed",
x="",
title = paste0("Ratio of Estimated Infections to Observed by State")) +
guides(color = guide_legend(override.aes = list(linewidth =2),
nrow=2,
byrow=TRUE))
ggsave(here('figures/ratio-estimated-to-observed-all-states.pdf'))
Comparing States at Specific Two-week Intervals
During Alpha Wave
plt1 <- subset %>%
# left_join(codes) %>%
filter(biweek == 7) %>%
ggplot(aes(y =fct_reorder(name,exp_cases_median/positive),
xmin = exp_cases_lb/positive,
xmax = exp_cases_ub/positive)) +
geom_errorbar() +
# geom_point(aes(y=state_name,x=exp_cases_median/positive),
# color = "darkred", size=.8, alpha=.5) +
theme_c() +
theme(axis.text.y=element_text(size=10),
axis.text.x = element_text(size=11),
axis.title = element_text(size = 17),
plot.subtitle = element_text(size=14),
plot.title = element_text(size=17, face="plain")) +
labs(y="", x="Ratio of Estimated Infections\nto Observed",
title ="During Two-week Period\nin the Alpha Wave",
subtitle = "March 26, 2021 through April 8, 2021")
plt1
During Delta Wave
subset %>%
mutate(ratio = exp_cases_median/positive) %>%
group_by(biweek) %>%
summarize(mean = mean(ratio),
mindate=min(date),
maxdate=max(date)) %>%
arrange(desc(mean))
plt2 <- subset %>%
# left_join(codes) %>%
filter(biweek == 13) %>%
ggplot(aes(y =fct_reorder(name,exp_cases_median/positive),
xmin = exp_cases_lb/positive,
xmax = exp_cases_ub/positive)) +
geom_errorbar() +
# geom_point(aes(y=state_name,x=exp_cases_median/positive),
# color = "darkred", size=.8, alpha=.5) +
theme_c() +
theme(axis.text.y=element_text(size=10),
axis.text.x = element_text(size=11),
axis.title = element_text(size = 17),
plot.subtitle = element_text(size=14),
plot.title = element_text(size=17, face="plain")) +
labs(y="", x="Ratio of Estimated Infections\nto Observed",
title ="During Two-week Period\nin the Delta Wave",
subtitle = "June 18, 2021 through July 1, 2021")
plt2
During Omicron Wave
subset %>%
left_join(codes) %>%
filter(biweek == 27) %>%
pull(date) %>%
range()
subset %>%
left_join(codes) %>%
filter(biweek == 13) %>%
pull(date) %>%
range()
plt3 <- subset %>%
# left_join(codes) %>%
filter(biweek == 27) %>%
ggplot(aes(y =fct_reorder(name,exp_cases_median/positive),
xmin = exp_cases_lb/positive,
xmax = exp_cases_ub/positive)) +
geom_errorbar() +
# geom_point(aes(y=state_name,x=exp_cases_median/positive),
# color = "darkred", size=.8, alpha=.5) +
theme_c() +
theme(axis.text.y=element_text(size=10),
axis.text.x = element_text(size=11),
axis.title = element_text(size = 17),
plot.subtitle = element_text(size=14),
plot.title = element_text(size=17, face="plain")) +
labs(y="", x="Ratio of Estimated Infections\nto Observed",
title ="During Two-week Period\nin the Omicron Wave",
subtitle = "December 31, 2021 through January 14, 2022")
plt3
All Waves Together
plt <- plot_grid(plt1,plt2, plt3, nrow=1, align="none")
title <- ggplot() +
labs(title = "Ratio of Estimated Infections to Observed Infections") +
theme_c(plot.title=element_text(size=25))
plot_grid(title, plt, nrow=2, rel_heights= c(.05,.95))
